package com.globopt.global.old;

public class Global {
	// input
	private int nvars;
	private Function fun;
	private double[] lb;
	private double[] ub;
	private Opts opts;

	// output
	private double[][] x0;
	private double[] f0;
	private int nc;
	private int nfe;

	public Global(Function fun) {
		this.fun = fun;
		this.nvars = fun.nvars;
		this.fun = new F1();
		this.opts = new Opts(nvars);

		this.lb = new double[nvars];
		this.ub = new double[nvars];
		for (int i = 0; i < nvars; i++) {
			this.lb[i] = -1.0;
			this.ub[i] = 1.0;
		}
	}

	public void calculate() throws Exception {
		// %======== Check number of best points selected ====%
		if (opts.getNg0() < 1)
			opts.setNg0(1);
		if (opts.getNg0() > 20)
			opts.setNg0(20);

		int ng10 = 100; // % maximum number of points selected
		int nnn = 0;
		int nn100;
		int n;

		// %======== Check number of sample points ========%

		if (opts.getN100() < 20) {
			opts.setN100(20);
		}
		if (opts.getN100() > 100000) {
			opts.setN100(100000);
		}
		if (opts.getN100() >= 1000) {
			nn100 = 1000;
			n = (int) Math.floor(opts.getN100() / 1000);
			nnn = opts.getN100() - n * 1000;
			opts.setN100(n * 1000);

		} else if (opts.getN100() < 1000 && opts.getN100() > 100) {
			nn100 = 100;
			n = (int) Math.floor(opts.getN100() / 100);
			nnn = opts.getN100() - n * 100;
			opts.setN100(n * 100);

		} else {
			nn100 = opts.getN100();
			n = 1;
		}
		// %======= Initialize the variables ========%

		// double[][] f = arraySetValues(100, 1, 0.0);
		this.f0 = new double[100];
		double[] f1 = new double[100];
		double[] fcl = new double[100];
		for (int i = 0; i < 100; i++) {
			this.f0[i] = 0.0;
			f1[i] = 0.0;
			fcl[i] = 0.0;
		}

		double[][] x = arraySetValues(nvars, 100, 0.0);
		this.x0 = arraySetValues(nvars, 100, 0.0);
		double[][] x1 = arraySetValues(nvars, 100, 0.0);
		double[][] xcl = arraySetValues(nvars, 100, 0.0);

		// double d = 2.0 * Math.sqrt(nvars) * 0.000001;
		double b1 = 1.0 / nvars;

		double[] f = new double[ng10];
		int[] ic = new int[ng10]; // % index of the cluster
		for (int i = 0; i < ng10; i++) {
			f[i] = 9.9e10;
			ic[i] = 0;
		}

		double[] pmax = new double[nvars];
		double[] pmin = new double[nvars];
		for (int i = 0; i < nvars; i++) {
			pmax[i] = (lb[i] - ub[i]) / 2.0;
			pmin[i] = lb[i] + ub[i];
		}

		// double alfa = 0.01;
		this.nfe = 0; // Number of function evaluations
		int ng = 0; // Current number of points selected
		int ns = 0; // Number of samplings
		this.nc = 0; // Number of local minimizers
		// int ncp = 1;
		int n0 = 0;
		int n1 = 0;
		int im = 0; // Index for FM
		double fm = 9.9e10;
		double relcon = Math.pow(10.0, -opts.getNsig());
		int it = 1;

		// %============ Main iteration =============%

		while (it > 0) {

			// %==================================================%
			// % SAMPLING %
			// %==================================================%

			n0 += opts.getN100();
			int nm = n0 - 1;
			ng += opts.getNg0();
			ns++;

			if (ns * opts.getNg0() > 100) {
				System.out.println("*** TOO MANY SAMPLING ***");
				printResults(pmin, pmax);
				return;
			}

			double b = Math.pow(1.0 - Math.pow(0.01, (1.0 / nm)), b1);
			double bb = 0.1 * b;

			double[] y = new double[this.nvars];
			for (int i = 0; i < n; i++) {
				double[][] r = MatlabOperationHelpers.twoDimensionalArraysFillWithRandomNumber(nn100, nvars);

				for (int j = 0; j < nn100; j++) {
					y = MatlabOperationHelpers.vectorAddition(MatlabOperationHelpers.vectorMultiply(2.0, r[j]), -1.0);
					double fc = MatlabOperationHelpers.globfun(y, fun, pmin, pmax);

					if (fc < fm) {
						f[im] = fc;
						for (int k = 0; k < this.nvars; k++) {
							x[k][im] = y[k];
						}
						ic[im] = 0;
						im = MatlabOperationHelpers.maxIndex(f);
						fm = f[im];
					}

				}
			}
			if (nnn > 0) {
				double[][] r = MatlabOperationHelpers.twoDimensionalArraysFillWithRandomNumber(nn100, nvars);
				for (int j = 0; j < nnn; j++) {
					y = MatlabOperationHelpers.vectorAddition(MatlabOperationHelpers.vectorMultiply(2.0, r[j]), -1.0);
					double fc = MatlabOperationHelpers.globfun(y, fun, pmin, pmax);
					if (fc < fm) {
						f[im] = fc;
						for (int k = 0; k < this.nvars; k++) {
							x[k][im] = y[k];
						}
						ic[im] = 0;
						im = MatlabOperationHelpers.maxIndex(f);
						fm = f[im];
					}
				}

			}

			nfe += opts.getN100() + nnn;

			if (opts.getDisplay().equals("on")) {
				System.out.println("Function evaluations used for sampling : " + opts.getN100());
			}

			// %==================================================%
			// % SORTING %
			// %==================================================%

			int[] is = MatlabOperationHelpers.sortIndex(f);
			f = MatlabOperationHelpers.sortByIndex(f, is);
			x = MatlabOperationHelpers.sortByIndex(x, is);
			ic = MatlabOperationHelpers.sortByIndex(ic, is);
			im = MatlabOperationHelpers.maxIndex(f);
			fm = f[im];

			// %==================================================%
			// % CLUSTERING TO X*
			// %==================================================%
			//
			// TODO
			if (nc > 0) {
				for (int iii = 0; iii < nc; iii++) {
					int in1 = 0;
					fcl[0] = f0[iii];
					xcl = MatlabOperationHelpers.copyColumn(xcl, 0, x0, iii);
					int[] icm = MatlabOperationHelpers.searchIndex(MatlabOperationHelpers.subArray(ic, 0, ng - 1), iii);
					if (icm.length > 0) {
						in1 += icm.length;
						for (int j = 1; j < in1 - 1; j++) {
							xcl = MatlabOperationHelpers.copyColumn(xcl, j, x, icm[j - 1]);
							fcl[j] = f[icm[j - 1]];
						}
					}
					slc(fcl, xcl, f, x, ic, ng, in1, 0, iii, b);
				}
			}

			// %==================================================%
			// % CLUSTERING TO X1 %
			// %==================================================%

			if (n1 > 0) {
				for (int iii = 0; iii < n1; iii++) {
					// int i = 0;
					// int in1 = i;
					fcl[0] = f1[iii];
					xcl = MatlabOperationHelpers.copyColumn(xcl, 0, x1, iii);
					slc(fcl, xcl, f, x, ic, ng, 0, 0, iii, b);
				}
			}

			// %==================================================%
			// % LOCAL SEARCH %
			// %==================================================%

			it = 0;
			int nsp = 0;
			double ff = Double.MIN_NORMAL;
			int nfevl = Integer.MIN_VALUE;
			double[] w = new double[nvars];
			int icj = Integer.MIN_VALUE;

			for (int i1 = 0; i1 < ng; i1++) {
				if (ic[i1] != 0) {
					continue;
				}
				for (int i = 0; i < nvars; i++) {
					y[i] = x[i][i1];
				}
				ff = f[i1];

				// %======== Call local solver =======%
				if (opts.getMethod().equals("bfgs")) {
					throw new Exception("bfgs not implemented in java");
				} else if (opts.getMethod().equals("unirandi")) {
					Unirandi unirandi = new Unirandi(fun, y, ff, pmin, pmax, opts.getMaxfn(), relcon);
					unirandi.calculate();
					y = unirandi.getX();
					ff = unirandi.getF();
					nfevl = unirandi.getNfev();

				} else {
					throw new Exception("wrong local search method");
				}

				if (nc > 0) {
					for (int iv = 0; iv < nc; iv++) {
						double a = 0.0;
						for (int i = 0; i < w.length; i++) {
							w[i] = Math.abs(x0[i][iv] - y[i]);
							if (w[i] > a) {
								a = w[i];
							}
						}
						if (a < bb) {
							n1++;
							if (opts.getDisplay().equals("on")) {
								System.out.println("New seed point added to the cluster no. " + iv + ", NFEV = " + nfevl);
								for (int i = 0; i < w.length; i++) {
									w[i] = x[i][i1] * pmax[i] + pmin[i];
								}
								System.out.println("ff= " + ff);
								System.out.println("w= " + MatlabOperationHelpers.arrayToString(w));
							}
							if (ff < f0[iv]) {
								if (opts.getDisplay().equals("on")) {
									System.out.println("*** IMPROVEMENT ON THE LOCAL MINIMUM NO. " + iv + " : " + f0[iv] + " FOR " + ff);
									w = MatlabOperationHelpers.vectorAddition(MatlabOperationHelpers.vectorMultiply(y, pmax), pmin);
									System.out.println("w= " + MatlabOperationHelpers.arrayToString(w));
								}
								f0[iv] = ff;
								for (int i = 0; i < nvars; i++) {
									x0[i][iv] = y[i];
								}
							}
							if (n1 > 20) {
								System.out.println("*** TWO MANY SEED POINTS ***");
								printResults(pmin, pmax);
								return;
							}
							x1 = MatlabOperationHelpers.copyColumn(x1, n1, x, i1);
							xcl = MatlabOperationHelpers.copyColumn(xcl, 0, x, i1);
							f1[n1] = f[i1];
							fcl[0] = f[i1];
							// ic1[n1] = iv; never used
							icj = iv;
							nsp = 1;
							break; // new seed point, go to clustering to the
									// new point
						}
					}
				}

				// %==================================================%
				// % NEW LOCAL MINIMUM %
				// %==================================================%

				if (nsp == 0) {

					nc++;
					// ncp++;
					if (opts.getDisplay().equals("on")) {
						System.out.println("'*** LOCAL MINIMUM NO. " + nc + " : " + ff + ", NFEV =" + nfevl);
						w = MatlabOperationHelpers.vectorAddition(MatlabOperationHelpers.vectorMultiply(y, pmax), pmin);
						System.out.println("W: " + MatlabOperationHelpers.arrayToString(w));
					}
					for (int i = 0; i < nvars; i++) {
						x0[i][nc] = y[i];
						xcl[i][0] = y[i];
					}
					fcl[0] = ff;
					f0[nc] = ff;
					if (nc >= opts.getMaxnlm()) {
						System.out.println("*** TOO MANY CLUSTERS ***");
						printResults(pmin, pmax);
						return;
					}
					it = 1;
					icj = nc;
				}

				// %==================================================%
				// % CLUSTERING TO THE NEW POINT %
				// %==================================================%

				nfe += nfevl;
				ic[i1] = icj;
				slc(fcl, xcl, f, x, ic, ng, 0, 0, icj, b);
				nsp = 0;
			}

		}
		if (opts.getDisplay().equals("on") || opts.getDisplay().equals("final")) {
			printResults(pmin, pmax);
		} else {
			sortF0X0();
			for (int i = 0; i < nc; i++) {
				for (int j = 0; j < nvars; j++) {
					x0[j][i] = x0[j][i] * pmax[j] + pmin[j];
				}
			}

		}
	}

	private double[][] arraySetValues(int row, int col, double value) {
		double[][] ret = new double[row][col];

		for (int i = 0; i < row; i++) {
			for (int j = 0; j < col; j++) {
				ret[i][j] = value;
			}
		}

		return ret;
	}

	private void sortF0X0() {
		for (int i = 0; i < this.nc; i++) {
			int minindex = i;
			for (int j = i + 1; j < this.nc; j++) {
				if (this.f0[j] < this.f0[minindex]) {
					minindex = j;
				}
			}
			double c;
			c = f0[i];
			f0[i] = f0[minindex];
			f0[minindex] = c;
			for (int j = 0; j < this.nvars; j++) {
				c = x0[j][i];
				x0[j][i] = x0[j][minindex];
				x0[j][minindex] = c;
			}
		}
	}

	private void printResults(double[] pmin, double[] pmax) {
		// %==================================================%
		// % PRINT RESULTS %
		// %==================================================%

		System.out.println(" NORMAL TERMINATION AFTER " + this.nfe + " FUNCTION EVALUATIONS");
		sortF0X0();
		for (int i = 0; i < this.nc; i++) {
			for (int j = 0; j < this.nvars; j++) {
				x0[j][i] = x0[j][i] * pmax[j] + pmin[j];

			}
		}

		System.out.println("LOCAL MINIMUM FOUND:" + nc);
		System.out.println("f0 = ");
		System.out.println(MatlabOperationHelpers.arrayToString(f0));
		System.out.println("x0 = ");
		System.out.println(MatlabOperationHelpers.arrayToString(x0));
		System.out.println("GLOBAL MINIMUM VALUE: " + f0[0]);
		System.out.println("GLOBAL MINIMUM: ");
		for (int i = 0; i < this.nvars; i++) {
			System.out.print(x0[i][0] + ", ");
		}
	}

	/**
	 * Single Linkage Clustering
	 */
	public void slc(double[] fcl, double[][] xcl, double[] f, double[][] x, int[] ic, int ng, int in1, int i, int iii, double b) {
		// %==================================================%
		// % Single Linkage Clustering %
		// %==================================================%

		while (i <= in1) {
			for (int j = 0; j < ng; j++) {

				if (ic[j] == 0 && fcl[i] < f[j]) {
					double[] w = new double[nvars];
					double a = 0.0;
					for (int k = 0; k < nvars; k++) {
						w[k] = Math.abs(xcl[k][i] - x[k][j]);
						if (w[k] > a) {
							a = w[k];
						}
					}
					if (a < b) {
						if (opts.getDisplay().equals("on")) {
							System.out.println("Sample point added to the cluster No. " + iii);
							double[] xc = new double[nvars];
							for (int k = 0; k < nvars; k++) {
								xc[k] = x[k][j] * ub[k] + lb[k];
							}
							System.out.println("F = ");
							System.out.println(MatlabOperationHelpers.arrayToString(f));
							System.out.println("X = ");
							System.out.println(MatlabOperationHelpers.arrayToString(xc));
						}
						if (in1 < 99) {
							in1++;
						}
						fcl[in1] = f[j];
						xcl = MatlabOperationHelpers.copyColumn(xcl, in1, x, j);
						ic[j] = iii;
					}
				}
			}
			i++;
		}
	}

	public Function getFun() {
		return fun;
	}

	public double[] getLb() {
		return lb;
	}

	public void setLb(double[] lb) {
		this.lb = lb;
	}

	public double[] getUb() {
		return ub;
	}

	public void setUb(double[] ub) {
		this.ub = ub;
	}

	public int getNvars() {
		return nvars;
	}

	public void setOpts(Opts opts) {
		this.opts = opts;
	}

	public Opts getOpts() {
		return opts;
	}

	public double[][] getX0() {
		return x0;
	}

	public double[] getF0() {
		return f0;
	}

	public int getNc() {
		return nc;
	}

	public int getNfe() {
		return nfe;
	}
}
